function [Beta,sigma2] = bayes_reg(Y,X,Beta,sigma2,prior,restr)
%Returns a draw from the posterior distribution of a linear regression
%   Detailed explanation goes here

      Treg = size(Y,1);
      K = size(X,2);
      
      % First get OLS Results
      
      Beta_OLS = inv(X'*X)*(X'*Y);
      SSE = (Y - X*Beta_OLS)'*(Y - X*Beta_OLS);
      sigma_OLS = SSE./(Treg-K);
      
      
      % Draw Beta conditional on Sigma
            
      if restr == 0
      
      % Set Priors
      
      if prior == 1

      Beta_0 = zeros(K,1);
      SigmaBeta_0 = 0.1*eye(K);
%     SigmaBeta_0(1,1) = 0.00000001;
            
      SigmaBeta_1 = inv(inv(SigmaBeta_0)+(inv(sigma2))*(X'*X));
      Beta_1 = SigmaBeta_1*(inv(SigmaBeta_0)*Beta_0+(inv(sigma2))*X'*Y);
      
      else if prior == 2
              
      Beta_0 = [zeros(1,K)]'; 
      
      tau = 0.01;
      decay = 2;    
    
      matrixOfPriorVars = tau./([1:K].^decay);
      SigmaBeta_0 = diag(matrixOfPriorVars);
      
      SigmaBeta_1 = inv(inv(SigmaBeta_0)+(inv(sigma2))*(X'*X));
      Beta_1 = SigmaBeta_1*(inv(SigmaBeta_0)*Beta_0+(inv(sigma2))*X'*Y);      
      
      else if prior == 0
              
      SigmaBeta_1 = sigma2*inv(X'*X);
      Beta_1 = Beta_OLS;
      
          end
          end
      end
      
      Beta = mvnrnd(Beta_1,SigmaBeta_1)';
      
      end
      
      % Draw Beta conditional on Sigma
      
       % Set Priors
      
      if prior == 1
          
      T0 = 1;
      D0 = 1;
      
      resids=Y-X*Beta;    % Compute residuals
      T1=T0+Treg;            %compute posterior df and scale matrix
      D1=D0+resids'*resids;
      z0=randn(T1,1);     %draw from IG
      z0z0=z0'*z0;
      sigma2=D1/z0z0;
      
      else if prior == 0
              
      S_post = SSE;
      v_post = Treg-K;
      sigma2 = inv(wish(inv(S_post),v_post));% Draw SIGMA
     
          end
      end
      
end

